Magnetostatics solution (nonlinear case)

%%capture
%run TEAM_13_current_J.ipynb
%run TEAM_13_nonlinearity.ipynb

##############################################################################
# Tree/Cotree gauging
##############################################################################
MESH = pde.mesh3.netgen(geoOCCmesh)
R = pde.tools.tree_cotree_gauge(MESH)

R

<219723x187864 sparse matrix of type '<class 'numpy.float64'>'
        with 187864 stored elements in Compressed Sparse Column format>

##############################################################################
# Assembly
##############################################################################

linear = '*coil,default'
nonlinear = 'r_steel,l_steel,mid_steel'
maxIter = 100

order = 1

fem_linear = pde.int.evaluate3(MESH, order = order, regions = linear).diagonal()
fem_nonlinear = pde.int.evaluate3(MESH, order = order, regions = nonlinear).diagonal()

D = pde.int.assemble3(MESH, order = order)

phix_Hcurl, phiy_Hcurl, phiz_Hcurl = pde.hcurl.assemble3(MESH, space = 'N0', matrix = 'M', order = order)
curlphix_Hcurl, curlphiy_Hcurl, curlphiz_Hcurl = pde.hcurl.assemble3(MESH, space = 'N0', matrix = 'K', order = order)

aJ = jx_L2 @ D @ phix_Hcurl.T +\
     jy_L2 @ D @ phiy_Hcurl.T +\
     jz_L2 @ D @ phiz_Hcurl.T

aJ = 1e7*aJ

def gss(A):
    curl_Ax = curlphix_Hcurl.T@A; curl_Ay = curlphiy_Hcurl.T@A; curl_Az = curlphiz_Hcurl.T@A

    Kxx = curlphix_Hcurl @ D @ sp.diags(fxx_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fxx_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear)@ curlphix_Hcurl.T
    Kyy = curlphiy_Hcurl @ D @ sp.diags(fyy_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fyy_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear)@ curlphiy_Hcurl.T
    Kzz = curlphiz_Hcurl @ D @ sp.diags(fzz_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fzz_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear)@ curlphiz_Hcurl.T

    Kxy = curlphiy_Hcurl @ D @ sp.diags(fxy_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fxy_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear)@ curlphix_Hcurl.T
    Kxz = curlphiz_Hcurl @ D @ sp.diags(fxz_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fxz_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear)@ curlphix_Hcurl.T

    Kyx = curlphix_Hcurl @ D @ sp.diags(fyx_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fyx_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear)@ curlphiy_Hcurl.T
    Kyz = curlphiz_Hcurl @ D @ sp.diags(fyz_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fyz_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear)@ curlphiy_Hcurl.T

    Kzx = curlphix_Hcurl @ D @ sp.diags(fzx_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fzx_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear)@ curlphiz_Hcurl.T
    Kzy = curlphiy_Hcurl @ D @ sp.diags(fzy_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fzy_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear)@ curlphiz_Hcurl.T

    return Kxx + Kyy + Kzz + Kxy + Kxz + Kyx + Kyz + Kzx + Kzy

def gs(A):
    curl_Ax = curlphix_Hcurl.T@A; curl_Ay = curlphiy_Hcurl.T@A; curl_Az = curlphiz_Hcurl.T@A
    return curlphix_Hcurl @ D @ (fx_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fx_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear) +\
           curlphiy_Hcurl @ D @ (fy_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fy_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear) +\
           curlphiz_Hcurl @ D @ (fz_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fz_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear) - aJ

def J(A):
    curl_Ax = curlphix_Hcurl.T@A; curl_Ay = curlphiy_Hcurl.T@A; curl_Az = curlphiz_Hcurl.T@A
    return np.ones(D.size)@ D @(f_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + f_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear) -aJ@A




import time
A = np.zeros(curlphix_Hcurl.shape[0])
mu = 0.0001
eps_newton = 1e-5
factor_residual = 1/2

tm = time.monotonic()
for i in range(maxIter):

    gssu = R.T @ gss(A) @ R
    gsu = R.T @ gs(A)

    tm = time.monotonic()
    wS = chol(gssu).solve_A(-gsu)
    # wS = sp.linalg.spsolve(gssu, -gsu)
    print('Solving took ', time.monotonic()-tm)

    w = R@wS

    alpha = 1

    # ResidualLineSearch
    # for k in range(1000):
    #     if np.linalg.norm(gs(A+alpha*w),np.inf) <= np.linalg.norm(gs(A),np.inf): break
    #     else: alpha = alpha*factor_residual

    # AmijoBacktracking
    float_eps = 1e-12; #float_eps = np.finfo(float).eps
    for kk in range(1000):
        if J(A+alpha*w)-J(A) <= alpha*mu*(gsu@wS) + np.abs(J(A))*float_eps: break
        else: alpha = alpha*factor_residual


    A_old_i = A
    A = A + alpha*w

    print ("NEWTON: Iteration: %2d " %(i+1)+"||obj: %2e" %J(A)+"|| ||grad||: %2e" %np.linalg.norm(R.T @ gs(A),np.inf)+"||alpha: %2e" % (alpha))


    # if ( np.linalg.norm(R.T @ gs(A),np.inf) < eps_newton):
    #     break
    if (np.abs(J(A)-J(A_old_i)) < 1e-5):
        break

elapsed = time.monotonic()-tm
print('Solving took ', elapsed, 'seconds')

C:\Users\Radu\AppData\Local\Temp\ipykernel_3464\2265775738.py:69: CholmodTypeConversionWarning:

converting matrix of class csr_matrix to CSC format

Solving took  9.32799999974668
NEWTON: Iteration:  1 ||obj: -2.103156e+10|| ||grad||: 1.787734e+07||alpha: 6.250000e-02
Solving took  9.921999999787658
NEWTON: Iteration:  2 ||obj: -1.705712e+11|| ||grad||: 1.735611e+07||alpha: 1.000000e+00
Solving took  9.858999999705702
NEWTON: Iteration:  3 ||obj: -1.814229e+11|| ||grad||: 1.481311e+07||alpha: 6.250000e-02
Solving took  9.844000000040978
NEWTON: Iteration:  4 ||obj: -1.876979e+11|| ||grad||: 1.235620e+07||alpha: 1.250000e-01
Solving took  9.734000000171363
NEWTON: Iteration:  5 ||obj: -1.893205e+11|| ||grad||: 1.173088e+07||alpha: 1.000000e+00
Solving took  10.125
NEWTON: Iteration:  6 ||obj: -1.914766e+11|| ||grad||: 2.433062e+07||alpha: 5.000000e-01
Solving took  10.125
NEWTON: Iteration:  7 ||obj: -1.926345e+11|| ||grad||: 1.811276e+07||alpha: 2.500000e-01
Solving took  9.75
NEWTON: Iteration:  8 ||obj: -1.932460e+11|| ||grad||: 1.037124e+07||alpha: 1.000000e+00
Solving took  10.015999999828637
NEWTON: Iteration:  9 ||obj: -1.935141e+11|| ||grad||: 1.238268e+07||alpha: 2.500000e-01
Solving took  9.891000000294298
NEWTON: Iteration: 10 ||obj: -1.937199e+11|| ||grad||: 8.582073e+06||alpha: 5.000000e-01
Solving took  10.030999999959022
NEWTON: Iteration: 11 ||obj: -1.937912e+11|| ||grad||: 6.299471e+06||alpha: 5.000000e-01
Solving took  9.875
NEWTON: Iteration: 12 ||obj: -1.938246e+11|| ||grad||: 6.385038e+06||alpha: 5.000000e-01
Solving took  9.875
NEWTON: Iteration: 13 ||obj: -1.938383e+11|| ||grad||: 6.572519e+06||alpha: 5.000000e-01
Solving took  9.875
NEWTON: Iteration: 14 ||obj: -1.938405e+11|| ||grad||: 6.023229e+06||alpha: 5.000000e-01
Solving took  9.82799999974668
NEWTON: Iteration: 15 ||obj: -1.938431e+11|| ||grad||: 6.974790e+06||alpha: 5.000000e-01
Solving took  9.82799999974668
NEWTON: Iteration: 16 ||obj: -1.938431e+11|| ||grad||: 6.077044e+06||alpha: 5.000000e-01
Solving took  10.32799999974668
NEWTON: Iteration: 17 ||obj: -1.938449e+11|| ||grad||: 6.955329e+06||alpha: 5.000000e-01
Solving took  10.485000000335276
NEWTON: Iteration: 18 ||obj: -1.938588e+11|| ||grad||: 9.856679e+05||alpha: 2.500000e-01
Solving took  10.155999999959022
NEWTON: Iteration: 19 ||obj: -1.938600e+11|| ||grad||: 2.636677e+05||alpha: 1.000000e+00
Solving took  10.046000000089407
NEWTON: Iteration: 20 ||obj: -1.938601e+11|| ||grad||: 3.536630e+04||alpha: 1.000000e+00
Solving took  10.030999999959022
NEWTON: Iteration: 21 ||obj: -1.938601e+11|| ||grad||: 1.316059e+03||alpha: 1.000000e+00
Solving took  10.094000000040978
NEWTON: Iteration: 22 ||obj: -1.938601e+11|| ||grad||: 9.937699e-01||alpha: 1.000000e+00
Solving took  10.063000000081956
NEWTON: Iteration: 23 ||obj: -1.938601e+11|| ||grad||: 2.234483e-03||alpha: 1.000000e+00
Solving took  10.032000000122935
NEWTON: Iteration: 24 ||obj: -1.938601e+11|| ||grad||: 4.892793e-06||alpha: 1.000000e+00
Solving took  10.04700000025332
NEWTON: Iteration: 25 ||obj: -1.938601e+11|| ||grad||: 2.767734e-06||alpha: 1.000000e+00
Solving took  10.109000000171363
NEWTON: Iteration: 26 ||obj: -1.938601e+11|| ||grad||: 2.256264e-06||alpha: 1.000000e+00
Solving took  10.094000000040978
NEWTON: Iteration: 27 ||obj: -1.938601e+11|| ||grad||: 2.483117e-06||alpha: 1.000000e+00
Solving took  10.20299999974668
NEWTON: Iteration: 28 ||obj: -1.938601e+11|| ||grad||: 3.411318e-06||alpha: 1.000000e+00
Solving took  10.203999999910593
NEWTON: Iteration: 29 ||obj: -1.938601e+11|| ||grad||: 2.939254e-06||alpha: 1.000000e+00
Solving took  10.483999999705702
NEWTON: Iteration: 30 ||obj: -1.938601e+11|| ||grad||: 2.426474e-06||alpha: 1.000000e+00
Solving took  11.061999999918044 seconds

curlphix_Hcurl_P0, curlphiy_Hcurl_P0, curlphiz_Hcurl_P0 = pde.hcurl.assemble3(MESH, space = 'N0', matrix = 'K', order = 0)
Bx = curlphix_Hcurl_P0.T @ A
By = curlphiy_Hcurl_P0.T @ A
Bz = curlphiz_Hcurl_P0.T @ A

##############################################################################
# Storing to vtk
##############################################################################

grid = pde.tools.vtklib.createVTK(MESH)
pde.tools.vtklib.add_L2_Vector(grid,Bx,By,Bz,'grad_x')
pde.tools.vtklib.writeVTK(grid, 'magnetostatics_solution.vtu')

import pyvista as pv
mesh = pv.read('magnetostatics_solution.vtu')
mesh2 = pv.read('current_density.vtu')

mesh.set_active_scalars("Scalars_")
threshed = mesh.threshold([0,4])

p = pv.Plotter()
p.add_mesh(threshed, style='surface', color = "w", opacity=0.2, label=None)

threshed.set_active_vectors("grad_x")
arrows = mesh.glyph(scale="grad_x", orient=True, tolerance=0.03, factor=18)
p.add_mesh(arrows, color="orange")


mesh2.set_active_vectors("J_L2")
arrows2 = mesh2.glyph(scale="J_L2", orient=True, tolerance=0.03, factor=9500.0)
p.add_mesh(arrows2, color="black")

p.camera_position = [(0, -500, 200),(0, 0, 0),(0, 0, 0)]
p.show(jupyter_backend="html")